#Replication code for "Accountability in Time: Evolution in Social-Accountability Institutions"
#by Brian Palmer-Rubin, Jésica E. Tapia Reyes, Daniel Berliner, Aaron Erlich, and Benjamin E. Bagozzi

# Analyses were performed using R version 4.3.2.

library(MASS)
library(lfe)
library(lmtest)
library(texreg)
library(data.table)
library(multiwayvcov)
library(interflex)
library(xtable)

setwd("YOUR FILE PATH HERE")


d1<-fread("WP_replication_final.csv")


# Description of variables included in the replication dataset:

# FOLIO: Unique ID for each request. Note that the request IDs match the FOLIO ID's used in the Infomex system, other than the year 2004, which did not have FOLIO IDs in Infomex’s public data. For the latter year, the FOLIO ID's correspond to the row number found in the original Infomex CSV file.
# date0: Request date
# YEAR: Year of request
# yearcount: Year of request, counting from 2003=0
# MONTH: Month of request
# month_ind: Combined month and year, e.g. July 2008 is 2008.5
# DEPENDENCIA: Dependencia (government agency to which the request was submitted)
# RESPUESTA: Official response designation
# MUNICIPIO: Municipality from which the request was filed
# CODIGOPOSTAL_NUMBER: Postcode from which the request was filed
# STATE: State from which the request was filed, three-letter abbreviation
# CODE_COMBINED: Combined municipality-state numeric code -- this is the unique identifier for municipalities
# MARGindex_2000: Index of economic marginality for each municipality, 2000 
# attach.inc: Indicator of requests with attached files included
# logwordcount: Logged count of words in the request text, including any attached files.
# readability: Modified readability score of the request text. Since many requests consist of only a single sentence, canonical readability metrics are not applicable here. Our readability score is thus calculated as the ratio of the number of characters to the number of words.
# medium_nonelec: An indicator for any requests filed non-electronically (and thus entered into the system by officials).
# LOGtimeWD1: Logged working days to response
# response_inexistencia: Responses that use the official designation claiming the requested information does not exist.
# S4_num_areas_prob: A measure of request complexity based on classifier predictions of whether the request had queries relating to multiple distinct administrative units of the agency.
# S4_num_info_queries_prob: A measure of request complexity based on classifier predictions of whether or not the request had multiple distinct queries.
# S5_distinct_reqs_related_prob: A measure of request complexity based on classifier predictions of whether the request contained distinct queries that were unrelated to each other. 
# S6_is_formal_prob: Classifier prediction of whether or not the request uses formal language
# S6_is_legal_prob: Classifier prediction of whether or not the request uses legal language
# S6_is_technical_prob: Classifier prediction of whether or not the request uses technical language
# S7_dummy_Data_prob: Classifier prediction of whether or not the request sought multiple data points.
# S7_dummy_Database_prob: Classifier prediction of whether or not the request sought a database. 
# S7_dummy_Datum_prob: Classifier prediction of whether or not the request sought a single data point.
# S7_dummy_Document_prob: Classifier prediction of whether or not the request sought a document.
# S7_dummy_MultipleDocuments_prob: Classifier prediction of whether or not the request sought multiple documents.
# S8_dummy_Activities_prob: Classifier prediction of whether or not the request pertained to institutional activities.
# S8_dummy_Budget_prob: Classifier prediction of whether or not the request pertained to budgets or spending.
# S8_dummy_Evaluation_prob: Classifier prediction of whether or not the request pertained to evaluations or statistics.
# S8_dummy_ExternalContracts_prob: Classifier prediction of whether or not the request pertained to external contracts.
# S8_dummy_InstStruc_prob: Classifier prediction of whether or not the request pertained to institutional structures, personnel, or human resources.
# S8_dummy_Other_prob: Classifier prediction of whether the request’s theme did not fall into the above themes.
# S8_dummy_Regulatory_prob: Classifier prediction of whether the request pertained to regulatory matters or permits.
# S10_is_clear_prob: Classifier prediction of whether the request is likely clear enough to respond to
# S10_is_competency_of_institution_prob: Classifier prediction of whether the request is likely to fall under the competency of the agency
# S10_is_public_prob: Classifier prediction of whether the request is likely to refer to information that is not classified
# S10_is_existant_prob: Classifier prediction of whether the request is likely to refer to information that exists in the records of the agency
# S11_dummy_Date_prob: Classifier prediction of whether the request mentions specific dates
# S11_dummy_Document_prob: Classifier prediction of whether the request mentions specific documents
# S11_dummy_Institution_prob: Classifier prediction of whether the request mentions specific institutions (inside of government)
# S11_dummy_Organization_prob: Classifier prediction of whether the request mentions specific organizations (outside of government)
# S11_dummy_Person_prob: Classifier prediction of whether the request mentions specific persons
# S11_dummy_Place_prob: Classifier prediction of whether the request mentions specific places
# Handcoded: Indicator for requests that were in the handcoded sample.
# R4_dummy_DeliveryElectronic: Classifier prediction of whether the response actually provides information in electronic form
# R4_dummy_Nonexistent: Classifier prediction of whether the response actually claims that the requested information does not exist
# quality: Classifier prediction of whether the response provides roughly half or more of the requested information
# HC_R_trueinex: For hand-coded sample only, indicator of responses that claimed the requested information did not exist. 
# HC_R_lessthanhalf: For hand-coded sample only, indicator of responses that provided less than half of the requested information.
# HC_R_noncomply: For hand-coded sample only, indicator of non-compliant responses. 
# HC_expert_combined: For hand-coded sample only, combined indicator of request expertise. 
# HC_inapproprequest: For hand-coded sample only, combined indicator of inappropriate requests. 
# HC_S4_num_info_queries: For hand-coded sample only, indicator of the same feature noted above.
# HC_S4_num_areas: For hand-coded sample only, indicator of the same feature noted above.
# HC_S5_distinct_reqs_related: For hand-coded sample only, indicator of the same feature noted above.
# HC_S7_dummy_Datum: For hand-coded sample only, indicator of the same feature noted above.
# HC_S7_dummy_Data: For hand-coded sample only, indicator of the same feature noted above.
# HC_S7_dummy_Database: For hand-coded sample only, indicator of the same feature noted above.
# HC_S7_dummy_Document: For hand-coded sample only, indicator of the same feature noted above.
# HC_S7_dummy_MultipleDocuments: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_Activities: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_Budget: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_Evaluation: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_ExternalContracts: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_InstStruc: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_Other: For hand-coded sample only, indicator of the same feature noted above.
# HC_S8_dummy_Regulator: For hand-coded sample only, indicator of the same feature noted above.

### Additional variables used in analyses but created in this replication code rather than included in the replication dataset:
# expert_combined: Sum of nine request measures; three are language-based and six are specificity-based. The language measures are classifier predictions of whether or not the request used formal, legal, or technical language. The specificity measures are classifier predictions of whether or not the request mentioned specific documents (by name or number), persons (by name or title), dates (more specific than a year), places (more specific than a state), governmental institutions (more specific than the name of the agency; for example a sub-secretariat or internal area), or organizations (private or non-governmental). The sum of the three language measures is divided by three, while the sum of the six language measures is divided by six, before the two sums are added together, so that the overall expertise score is comprised equally of language and specificity measures.
# expert_PCA:  As a robustness check, we use the first dimension of the results of a principal components analysis of the nine language and specificity measures discussed for the measure "expert_combined."
# R_trueinex: Indicator of whether or not the response actually claims the information requested does not exist, regardless of whether that denial is made through the official designation or in the text of the response letter (based on classifier prediction).
# R_lessthanhalf: Indicator of whether or not the response provides less than half of the requested information, based on combination of officially-designated negative response categories and classifier predictions of the response texts for officially-designated positive responses. 
# R_noncomply: Indicator of whether or not the response is legally noncompliant, based on comparison of the official response designation and classifier predictions of whether or not information was actually provided. 
# agency_expert_combined_mean_lag: Average expertise score of all requests received by a given agency in prior year.
# agency_expertPCA_mean_lag: Average expertise score of all requests received by a given agency in prior year (using alternative PCA measure of expertise). 
# top100DEPin_any_year: Indicator of agencies that were in the top 100 by request volume for any given year in the sample.



### Making expertise measure and other independent variables using classifier-predicted measures

d1$expert<-d1$S6_is_formal_prob + d1$S6_is_legal_prob + d1$S6_is_technical_prob
d1$specific<-d1$S11_dummy_Date_prob + d1$S11_dummy_Document_prob + d1$S11_dummy_Institution_prob + d1$S11_dummy_Organization_prob + d1$S11_dummy_Person_prob + d1$S11_dummy_Place_prob
d1$expert_combined<- (d1$expert/3 + d1$specific/6)

pca_output<-prcomp(d1[,c("S6_is_formal_prob","S6_is_legal_prob","S6_is_technical_prob","S11_dummy_Date_prob", "S11_dummy_Document_prob","S11_dummy_Institution_prob","S11_dummy_Organization_prob","S11_dummy_Person_prob","S11_dummy_Place_prob")])
cor(d1$expert_combined, pca_output$x[,1])
#orient the PCA measure in the same direction as the additive one
d1$expert_PCA<- (-1)*pca_output$x[,1]

# measure of "appropriateness" of the request e.g. clarity, sent to the right agency, info likely to exist and to not be classified
# inappropriate request = sum of all where those are lacking
d1$inapproprequest<- (1-d1$S10_is_clear_prob) + (1-d1$S10_is_competency_of_institution_prob) + (1-d1$S10_is_public_prob) + (1-d1$S10_is_existant_prob)




### Making response measures using classifier predictions

# True inexistencia: Response claims that the information does not exist, combining official designations with predictions based on response texts
# NAs are those where information was delivered electronically but do not have any value for the predictions
d1$R_trueinex<-NA
d1$R_trueinex[d1$RESPUESTA!="Inexistencia de la información solicitada" & d1$RESPUESTA!="Entrega de información en medio electrónico"]<-0
d1$R_trueinex[d1$RESPUESTA=="Inexistencia de la información solicitada"]<-1
d1$R_trueinex[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_Nonexistent>=0.5]<-1
d1$R_trueinex[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_Nonexistent<0.5]<-0

# Uninformative response: Either information not delivered or less than half of the requested information was provided, combining official designations with predictions based on response texts
# NAs are those where information was delivered electronically but do not have any value for the predictions
d1$R_infoprovided<-NA
d1$R_infoprovided[d1$RESPUESTA!="Entrega de información en medio electrónico"]<-0
d1$R_infoprovided[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_DeliveryElectronic>=0.5 & d1$quality>=0.5]<-1
d1$R_infoprovided[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_DeliveryElectronic<0.5]<-0
d1$R_infoprovided[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$quality<0.5]<-0
d1$R_lessthanhalf<- 1 - d1$R_infoprovided

# Noncompliant response: Where official designation claims information delivered but classifier prediction is that this is not the case.
# Here also NA anywhere that the official designation was not "Entrega", because legal noncompliance is only in relation to that.
d1$R_noncomply<-NA
d1$R_noncomply[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_DeliveryElectronic<0.5]<-1
d1$R_noncomply[d1$RESPUESTA=="Entrega de información en medio electrónico" & d1$R4_dummy_DeliveryElectronic>=0.5]<-0




### Make lagged agency expertise score measures
### NOTE: This can take roughly 20 minutes to run. 

d1_temp<-d1[,c("FOLIO","YEAR","DEPENDENCIA","expert_combined","expert_PCA")]

## create an agency x year dataframe
ud<-unique(d1_temp$DEPENDENCIA)
df1<- as.data.frame(cbind(
  rep(ud, each=length(2003:2019)),
  rep(2003:2019, length(ud))
  ))
names(df1)<-c("DEPENDENCIA","YEAR")
df1$YEAR<-as.numeric(df1$YEAR)


df1$expert_combined_mean<-NA
df1$expertPCA_mean<-NA


timestamp()
for(i in 1:nrow(df1)){
  tempdata<-d1_temp[d1_temp$DEPENDENCIA==df1$DEPENDENCIA[i] & d1_temp$YEAR==df1$YEAR[i],]
  df1$expert_combined_mean[i]<-mean(tempdata$expert_combined)
  df1$expertPCA_mean[i]<-mean(tempdata$expert_PCA)

  if(i %in% seq(100,14000,100)) print(i)
}
timestamp()


#lagging variables

lagger<-function(variable, country, year, laglength){
  
  country<-as.character(country)
  laggedvar<-rep(NA,length(variable))
  
  leadingNAs<-rep(NA,laglength)
  countryshift<-c(leadingNAs, country[1:(length(country)-laglength)])
  
  variableshift<-c(leadingNAs, variable[1:(length(variable)-laglength)])
  
  replacementrefs<-country==countryshift
  replacementrefs[is.na(replacementrefs)==T]<-FALSE
  laggedvar[replacementrefs]<-variableshift[replacementrefs]
  
  laggedvar 
  }

df1$expert_combined_mean_lag<-lagger(df1$expert_combined_mean, df1$DEPENDENCIA, df1$YEAR, 1)
df1$expertPCA_mean_lag<-lagger(df1$expertPCA_mean, df1$DEPENDENCIA, df1$YEAR, 1)

d1_temp$agency_expert_combined_mean_lag<-NA
d1_temp$agency_expertPCA_mean_lag<-NA

# inserting new variable into main dataset
timestamp()
for(i in 1:nrow(df1)){
  if(!is.na(df1$expert_combined_mean_lag[i])) d1_temp$agency_expert_combined_mean_lag[d1_temp$YEAR== df1$YEAR[i] & d1_temp$DEPENDENCIA==df1$DEPENDENCIA[i]] <- df1$expert_combined_mean_lag[i]
  if(!is.na(df1$expertPCA_mean_lag[i])) d1_temp$agency_expertPCA_mean_lag[d1_temp$YEAR== df1$YEAR[i] & d1_temp$DEPENDENCIA==df1$DEPENDENCIA[i]] <- df1$expertPCA_mean_lag[i]
  if(i %in% seq(100,14000,100)) print(i)
}
timestamp()

summary(d1_temp$agency_expert_combined_mean_lag)
summary(d1_temp$agency_expertPCA_mean_lag)

dim(d1)
d1<-merge(d1, d1_temp[,c("FOLIO","agency_expert_combined_mean_lag","agency_expertPCA_mean_lag")], by="FOLIO")
dim(d1)



### Indicator for main sample inclusion: all agencies that were among the top 100-requested agencies for at least one year in our sample period. 

top100in_any_year<-NULL

for(i in 2003:2019){
  top100in_any_year<-c(top100in_any_year, names(sort(table(d1$DEPENDENCIA[d1$YEAR==i]),dec=T)[1:100]))
}

length(unique(top100in_any_year))
top100DEPin_any_year<-unique(top100in_any_year)
d1$top100DEPin_any_year<-as.numeric(d1$DEPENDENCIA %in% top100DEPin_any_year)





### Descriptive figures and summary statistics table


# Figure 1: Trends over time in overall volume of information requests (panel 1) and in the average expertise score of information requests (panel 2), both by month.

pdf(file="Figure1.pdf", height=12, width=9)

par(mfrow=c(2,1))

month_temp<-tapply(d1$month_ind[d1$YEAR<2020], d1$month_ind[d1$YEAR<2020], length)

plot(as.numeric(names(month_temp)),month_temp, ylim=c(0, 26000),type="n",
   xlab="",ylab="Requests per Month", main="Monthly Request Volume Over Time", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)

l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)

axis(side=1, at=2003:2019)

legend(2007.6, 27100, lty=c(1,2), lwd=c(2,1), col=c("gray70","black"), legend=c("New President in Office","Constitutional/Legal Change"),
  bty="n",cex=0.8)


month_temp<-tapply(d1$expert_combined[d1$YEAR<2020], d1$month_ind[d1$YEAR<2020], mean)

plot(as.numeric(names(month_temp)),month_temp, ylim=c(0.145, 0.3),type="n",
   xlab="",ylab="Average Expertise Score", main="Monthly Average Expertise Score Over Time", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)

l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)

axis(side=1, at=2003:2019)

legend(2007.6, 0.305, lty=c(1,2), lwd=c(2,1), col=c("gray70","black"), legend=c("New President in Office","Constitutional/Legal Change"),
  bty="n",cex=0.8)

dev.off()



# descriptive statistics on magnitude of changes in expert score
(mean(d1$expert_combined[d1$YEAR==2019]) - mean(d1$expert_combined[d1$YEAR==2003])) / sd(d1$expert_combined)
mean( d1$expert_combined[d1$YEAR==2003]<= mean(d1$expert_combined[d1$YEAR==2019]) )





#Figure C.2 Histograms of the relative distributions of expertise score in 2003 and 2019

hist2003<-NULL
hist2019<-NULL

hist2003 <- hist(d1$expert_combined[d1$YEAR==2003], breaks=40, plot = FALSE)
hist2003$density <- hist2003$counts/sum(hist2003$counts)*100

hist2019 <- hist(d1$expert_combined[d1$YEAR==2019], breaks=40, plot = FALSE)
hist2019$density <- hist2019$counts/sum(hist2019$counts)*100

pdf("FigureC2.pdf", height=7, width=10)
plot(hist2003, freq=FALSE, xlim=c(0,1.1), ylim=c(0,16), col="gray60", xlab="Expertise Score", 
     ylab="Density", main="Distributions of Expertise Score in 2003 and 2019")
plot(hist2019, freq=FALSE, xlim=c(0,1.1), ylim=c(0,16), density=30, angle=60, col="black", add=T)
legend("topright", legend=c("2003","2019"), col=c("gray60","black"), density=c(NA,30), angle=c(NA,60))
dev.off()





# Summary statistics table

d1$year_ind<-d1$month_ind - 2003

variablelist<-c(
  "LOGtimeWD1","response_inexistencia","R_trueinex","R_lessthanhalf","R_noncomply",
  "expert_combined","year_ind","agency_expert_combined_mean_lag",
  "logwordcount","readability","attach.inc","medium_nonelec","inapproprequest",
  "S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob",
  "S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob",
  "S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Regulatory_prob","S8_dummy_Other_prob"
  )

sumstats<-round( cbind(
  apply(d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1, ..variablelist],2,min,na.rm=T),
  apply(d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1, ..variablelist],2,max,na.rm=T),
  apply(d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1, ..variablelist],2,mean,na.rm=T),
  apply(d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1, ..variablelist],2,median,na.rm=T),
  apply(d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1, ..variablelist],2,sd,na.rm=T) ), 3)

colnames(sumstats)<-c("Min.","Max.","Mean","Median","SD")
xtable(sumstats)




# Figure D.3: Trends over time in five different measures of responsiveness to information requests, by monthly averages.

pdf("FigureD3.pdf", height=8, width=13)

par(mfrow=c(2,3))

month_temp<-tapply(d1$LOGtimeWD1[d1$RESPUESTA!="Sin respuesta"], d1$month_ind[d1$RESPUESTA!="Sin respuesta"], mean, na.rm=T)
plot(as.numeric(names(month_temp)),month_temp, #ylim=c(0.145, 0.3),
  type="n", xlab="",ylab="Avg. Log Days-to-Response", main="Monthly Average Log Days-to-Response", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)
l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)
axis(side=1, at=2003:2019)

month_temp<-tapply(d1$response_inexistencia[d1$RESPUESTA!="Sin respuesta"], d1$month_ind[d1$RESPUESTA!="Sin respuesta"], mean, na.rm=T)
plot(as.numeric(names(month_temp)),month_temp, #ylim=c(0.145, 0.3),
  type="n", xlab="",ylab="Prop. Official Inex.", main="Monthly Proportion Official Inexistencia", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)
l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)
axis(side=1, at=2003:2019)

month_temp<-tapply(d1$R_trueinex[d1$RESPUESTA!="Sin respuesta"], d1$month_ind[d1$RESPUESTA!="Sin respuesta"], mean, na.rm=T)
plot(as.numeric(names(month_temp)),month_temp, #ylim=c(0.145, 0.3),
  type="n", xlab="",ylab="Prop. True Inex.", main="Monthly Proportion True Inexistencia", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)
l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)
axis(side=1, at=2003:2019)

month_temp<-tapply(d1$R_lessthanhalf[d1$RESPUESTA!="Sin respuesta"], d1$month_ind[d1$RESPUESTA!="Sin respuesta"], mean, na.rm=T)
plot(as.numeric(names(month_temp)),month_temp, #ylim=c(0.145, 0.3),
  type="n", xlab="",ylab="Prop. Uninformative", main="Monthly Proportion Uninformative Response", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)
l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)
axis(side=1, at=2003:2019)

month_temp<-tapply(d1$R_noncomply[d1$RESPUESTA!="Sin respuesta"], d1$month_ind[d1$RESPUESTA!="Sin respuesta"], mean, na.rm=T)
plot(as.numeric(names(month_temp)),month_temp, #ylim=c(0.145, 0.3),
  type="n", xlab="",ylab="Prop. Noncompliant", main="Monthly Proportion Noncompliant Response", xaxt="n")
abline(v=c(2006.91666666667, 2012.91666666667, 2018.91666666667), lwd=3, col="gray70")
abline(v=c(2007.5, 2014.08333333333, 2015.33333333333), lwd=1, lty=2)
lines(as.numeric(names(month_temp)),month_temp, lwd=1)
l1 <- loess(month_temp ~ as.numeric(names(month_temp)))
lines(as.numeric(names(month_temp)),l1$fitted,lwd=2)
axis(side=1, at=2003:2019)

dev.off()






### Analyses of request expertise


#Figure B.1: Comparisons of Request Expertise Score with municipality-level Economic Marginality

munis_marg<-tapply(d1$MARGindex_2000[d1$YEAR<2020], d1$MUNICIPIO[d1$YEAR<2020], mean, na.rm=T)
munis_expert<-tapply(d1$expert_combined[d1$YEAR<2020], d1$MUNICIPIO[d1$YEAR<2020], mean, na.rm=T)
munis_num<-tapply(d1$expert_combined[d1$YEAR<2020], d1$MUNICIPIO[d1$YEAR<2020], length)

dels_expert<-tapply(d1$expert_combined[d1$YEAR<2020 & d1$STATE=="DIF" & d1$MUNICIPIO!="TUXPAN"], d1$MUNICIPIO[d1$YEAR<2020 & d1$STATE=="DIF" & d1$MUNICIPIO!="TUXPAN"], mean)
dels_marg<-tapply(d1$MARGindex_2000[d1$YEAR<2020 & d1$STATE=="DIF" & d1$MUNICIPIO!="TUXPAN"], d1$MUNICIPIO[d1$YEAR<2020 & d1$STATE=="DIF" & d1$MUNICIPIO!="TUXPAN"], mean)
# Note there there are a lone 5 requests from 2008-2009 appearing as filed from "Tuxpan" and the Federal District, even though Tuxpan is not one of the delegaciones of the Federal District. These are removed here so as to not result in an extra point in the scatterplot.
dels_names<-names(dels_expert)
dels_names[dels_names=="MAGDALENA CONTRERAS  LA"]<-"LA MAGDALENA CONTRERAS"


pdf("FigureB1.pdf", height=12, width=9)

par(mfrow=c(2,1))

plot(log(munis_marg), munis_expert, xlab="Log Economic Marginality", ylab="Average Request Expertise", type="n",
  main="Local Econ. Marginality and Request Expertise by Municipality")
points(log(munis_marg), munis_expert, cex=log(munis_num+1)/3)
abline(lm(munis_expert~log(munis_marg)), lty=3)

plot(log(dels_marg), dels_expert, xlab="Log Economic Marginality", ylab="Average Request Expertise", type="n", xlim=c(1.5,3),
  main="Local Econ. Marginality and Request Expertise for Delegations within the Federal District")
text(log(dels_marg), dels_expert, labels=dels_names, cex=0.7)
abline(lm(dels_expert~log(dels_marg)), lty=3)

dev.off()




#Table 1: Linear models of request expertise, for requests 2003-2016.

d1$year_ind<-d1$month_ind - 2003
d1$maxTheme<-paste("Theme",as.character(apply(d1[,c("S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Regulatory_prob")],1,which.max)),sep="_")

# increasingly fine-grained fixed effects
d1$munipost<-paste(d1$CODE_COMBINED, d1$CODIGOPOSTAL_NUMBER, sep="_")
d1$munidep<-paste(d1$CODE_COMBINED, d1$DEPENDENCIA, sep="_")
d1$munipostdep<-paste(d1$CODE_COMBINED, d1$CODIGOPOSTAL_NUMBER, d1$DEPENDENCIA, sep="_")
d1$munideptheme<-paste(d1$CODE_COMBINED, d1$DEPENDENCIA,  d1$maxTheme, sep="_")
d1$munipostdeptheme<-paste(d1$CODE_COMBINED, d1$CODIGOPOSTAL_NUMBER, d1$DEPENDENCIA, d1$maxTheme, sep="_")


m1<-felm(expert_combined ~ year_ind | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m2<-felm(expert_combined ~ year_ind | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m3<-felm(expert_combined ~ year_ind | munipost | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m4<-felm(expert_combined ~ year_ind | munidep | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m5<-felm(expert_combined ~ year_ind | munipostdep | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m6<-felm(expert_combined ~ year_ind | munideptheme | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])
m7<-felm(expert_combined ~ year_ind | munipostdeptheme | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d1[d1$YEAR<2017 & !is.na(d1$CODE_COMBINED) & !is.na(d1$CODIGOPOSTAL_NUMBER),])

texreg(list(m1, m2, m3, m4, m5, m6, m7), digits=4,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))




### Analyses of responsiveness

d2<-d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1,]


# Table 2

m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    



# Table 3

m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    








### Robustness checks

#Tables I.1 and I.2: Robustness checks with alternate PCA version of Expert Request Score

m1 <- felm(LOGtimeWD1 ~ 
  expert_PCA*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_PCA*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_PCA*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_PCA*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_PCA*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    




m1 <- felm(LOGtimeWD1 ~ 
  expert_PCA*agency_expertPCA_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_PCA*agency_expertPCA_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_PCA*agency_expertPCA_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_PCA*agency_expertPCA_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_PCA*agency_expertPCA_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d2, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    






#Tables I.3 and I.4: Robustness check with larger sample including all agencies that appear in the data.

d3<-d1[d1$RESPUESTA!="Sin respuesta",]

m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    


m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d3, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))    

rm(d3)



# Tables I.7 and I.8: Robustness checks including additional controls for time-periods

d1$sexenio1<-"Fox"
d1$sexenio1[d1$month_ind>=2006.91666666667]<-"Calderon"
d1$sexenio1[d1$month_ind>=2012.91666666667]<-"EPN"
d1$sexenio1[d1$month_ind>=2018.91666666667]<-"AMLO"

d1$policy_periods<-1
d1$policy_periods[d1$month_ind>=2007.5]<-2
d1$policy_periods[d1$month_ind>=2014.08333333333]<-3
d1$policy_periods[d1$month_ind>=2015.33333333333]<-4

table(d1$sexenio1)
table(d1$policy_periods)

d4<-d1[d1$RESPUESTA!="Sin respuesta" & d1$top100DEPin_any_year==1,]

m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*year_ind+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*year_ind+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*year_ind+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*year_ind+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*year_ind+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1))    



m1 <- felm(LOGtimeWD1 ~ 
  expert_combined*agency_expert_combined_mean_lag+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m2 <- felm(response_inexistencia ~ 
  expert_combined*agency_expert_combined_mean_lag+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m3 <- felm(R_trueinex ~ 
  expert_combined*agency_expert_combined_mean_lag+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m4 <- felm(R_lessthanhalf ~ 
  expert_combined*agency_expert_combined_mean_lag+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

m5 <- felm(R_noncomply ~ 
  expert_combined*agency_expert_combined_mean_lag+as.factor(sexenio1)+as.factor(policy_periods)+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) | DEPENDENCIA | 0 | DEPENDENCIA, data=d4, exactDOF=TRUE,psdef=FALSE)

texreg(list(m1, m2, m3, m4, m5), digits=3,  stars = c(0.01, 0.05, 0.1))    

rm(d4)



#Tables I.5 and I.6: Robustness checks using hand-coded sample only.

m1 <- lm(LOGtimeWD1 ~ 
  HC_expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m1_se<- cluster.vcov(m1, d1$DEPENDENCIA)

m2 <- lm(HC_R_trueinex ~ 
  HC_expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m2_se<- cluster.vcov(m2, d1$DEPENDENCIA)

m3 <- lm(HC_R_lessthanhalf ~ 
  HC_expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m3_se<- cluster.vcov(m3, d1$DEPENDENCIA)

m4 <- lm(HC_R_noncomply ~ 
  HC_expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m4_se<- cluster.vcov(m4, d1$DEPENDENCIA)

texreg(list(m1, m2, m3, m4), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"), 
  override.se=list(coeftest(m1, vcov.=m1_se)[,2],
                   coeftest(m2, vcov.=m2_se)[,2],
                   coeftest(m3, vcov.=m3_se)[,2],
                   coeftest(m4, vcov.=m4_se)[,2]), 
  override.p=list(coeftest(m1, vcov.=m1_se)[,4],
                  coeftest(m2, vcov.=m2_se)[,4],
                  coeftest(m3, vcov.=m3_se)[,4],
                  coeftest(m4, vcov.=m4_se)[,4]))    


m1 <- lm(LOGtimeWD1 ~ 
  HC_expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m1_se<- cluster.vcov(m1, d1$DEPENDENCIA)

m2 <- lm(HC_R_trueinex ~ 
  HC_expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m2_se<- cluster.vcov(m2, d1$DEPENDENCIA)

m3 <- lm(HC_R_lessthanhalf ~ 
  HC_expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m3_se<- cluster.vcov(m3, d1$DEPENDENCIA)

m4 <- lm(HC_R_noncomply ~ 
  HC_expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec +HC_inapproprequest + HC_S4_num_info_queries + HC_S4_num_areas + HC_S5_distinct_reqs_related +
  HC_S7_dummy_Datum + HC_S7_dummy_Data + HC_S7_dummy_Database + HC_S7_dummy_Document + HC_S7_dummy_MultipleDocuments + HC_S8_dummy_Activities + HC_S8_dummy_Budget + HC_S8_dummy_Evaluation + HC_S8_dummy_ExternalContracts + HC_S8_dummy_InstStruc + HC_S8_dummy_Other + HC_S8_dummy_Regulatory, 
  data=d1)
m4_se<- cluster.vcov(m4, d1$DEPENDENCIA)

texreg(list(m1, m2, m3, m4), digits=3,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"), 
  override.se=list(coeftest(m1, vcov.=m1_se)[,2],
                   coeftest(m2, vcov.=m2_se)[,2],
                   coeftest(m3, vcov.=m3_se)[,2],
                   coeftest(m4, vcov.=m4_se)[,2]), 
  override.p=list(coeftest(m1, vcov.=m1_se)[,4],
                  coeftest(m2, vcov.=m2_se)[,4],
                  coeftest(m3, vcov.=m3_se)[,4],
                  coeftest(m4, vcov.=m4_se)[,4]))    








### Visualizations of primary model results from Tables 2 and 3

# Repeating the models from Tables 2 and 3 but using lm() rather than felm() for plotting purposes. Note that numeric results are identical. 

A1 <- lm(LOGtimeWD1 ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
A1_se<- cluster.vcov(A1, d2$DEPENDENCIA)

A2 <- lm(response_inexistencia ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
A2_se<- cluster.vcov(A2, d2$DEPENDENCIA)

A3 <- lm(R_trueinex ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
A3_se<- cluster.vcov(A3, d2$DEPENDENCIA)

A4 <- lm(R_lessthanhalf ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
A4_se<- cluster.vcov(A4, d2$DEPENDENCIA)

A5 <- lm(R_noncomply ~ 
  expert_combined*year_ind+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
A5_se<- cluster.vcov(A5, d2$DEPENDENCIA)




expert_hi <- mean(d2$expert_combined) + sd(d2$expert_combined)
expert_lo <- mean(d2$expert_combined) - sd(d2$expert_combined)

# confirm number of group fixed effects
length(names(coef(A1))[grep("DEPENDENCIA",names(coef(A1)))]) #208
length(names(coef(A2))[grep("DEPENDENCIA",names(coef(A2)))]) #208
length(names(coef(A4))[grep("DEPENDENCIA",names(coef(A4)))]) #208
length(names(coef(A5))[grep("DEPENDENCIA",names(coef(A5)))]) #208


predfunction_year<-function(model, vcov){

    pe<-coef(model)
    vc<-vcov
    #vc<-cluster.vcov(model, d2$DEPENDENCIA)
    sims <- 1000
    simbetas <- mvrnorm(sims,pe,vc)

    valueseq<- seq(0.5,16.5,1)
    EVs_LO<-rep(NA,length(valueseq))
    lower_LO<-rep(NA,length(valueseq))
    upper_LO<-rep(NA,length(valueseq))
    EVs_HI<-rep(NA,length(valueseq))
    lower_HI<-rep(NA,length(valueseq))
    upper_HI<-rep(NA,length(valueseq))

    for (i in 1:length(valueseq)){
        xhyp_LO<-c( 1, 
        expert_combined=expert_lo,
        year_ind=valueseq[i],
        logwordcount = mean(d2$logwordcount),
        readability = mean(d2$readability),
        attach.inc = mean(d2$attach.inc),
        medium_nonelec = mean(d2$medium_nonelec),
        inapproprequest = mean(d2$inapproprequest),
        S4_num_info_queries_prob = mean(d2$S4_num_info_queries_prob),
        S4_num_areas_prob = mean(d2$S4_num_areas_prob),
        S5_distinct_reqs_related_prob = mean(d2$S5_distinct_reqs_related_prob),
        S7_dummy_Datum_prob = mean(d2$S7_dummy_Datum_prob),
        S7_dummy_Data_prob = mean(d2$S7_dummy_Data_prob),
        S7_dummy_Database_prob = mean(d2$S7_dummy_Database_prob),
        S7_dummy_Document_prob = mean(d2$S7_dummy_Document_prob),
        S7_dummy_MultipleDocuments_prob = mean(d2$S7_dummy_MultipleDocuments_prob),
        S8_dummy_Activities_prob = mean(d2$S8_dummy_Activities_prob),
        S8_dummy_Budget_prob = mean(d2$S8_dummy_Budget_prob),
        S8_dummy_Evaluation_prob = mean(d2$S8_dummy_Evaluation_prob),
        S8_dummy_ExternalContracts_prob = mean(d2$S8_dummy_ExternalContracts_prob),
        S8_dummy_InstStruc_prob = mean(d2$S8_dummy_InstStruc_prob),
        S8_dummy_Other_prob = mean(d2$S8_dummy_Other_prob),
        S8_dummy_Regulatory_prob = mean(d2$S8_dummy_Regulatory_prob),  
        rep(1/11, 11),
        rep(1/208, 208),
        valueseq[i] * expert_lo
        )
      
      pred<- simbetas %*% xhyp_LO
      EVs_LO[i]<-quantile(pred,0.5)
      lower_LO[i]<-quantile(pred,0.025)
      upper_LO[i]<-quantile(pred,0.975)   } 


    for (i in 1:length(valueseq)){
        xhyp_HI<-c( 1, 
        expert_combined=expert_hi,
        year_ind=valueseq[i],
        logwordcount = mean(d2$logwordcount),
        readability = mean(d2$readability),
        attach.inc = mean(d2$attach.inc),
        medium_nonelec = mean(d2$medium_nonelec),
        inapproprequest = mean(d2$inapproprequest),
        S4_num_info_queries_prob = mean(d2$S4_num_info_queries_prob),
        S4_num_areas_prob = mean(d2$S4_num_areas_prob),
        S5_distinct_reqs_related_prob = mean(d2$S5_distinct_reqs_related_prob),
        S7_dummy_Datum_prob = mean(d2$S7_dummy_Datum_prob),
        S7_dummy_Data_prob = mean(d2$S7_dummy_Data_prob),
        S7_dummy_Database_prob = mean(d2$S7_dummy_Database_prob),
        S7_dummy_Document_prob = mean(d2$S7_dummy_Document_prob),
        S7_dummy_MultipleDocuments_prob = mean(d2$S7_dummy_MultipleDocuments_prob),
        S8_dummy_Activities_prob = mean(d2$S8_dummy_Activities_prob),
        S8_dummy_Budget_prob = mean(d2$S8_dummy_Budget_prob),
        S8_dummy_Evaluation_prob = mean(d2$S8_dummy_Evaluation_prob),
        S8_dummy_ExternalContracts_prob = mean(d2$S8_dummy_ExternalContracts_prob),
        S8_dummy_InstStruc_prob = mean(d2$S8_dummy_InstStruc_prob),
        S8_dummy_Other_prob = mean(d2$S8_dummy_Other_prob),
        S8_dummy_Regulatory_prob = mean(d2$S8_dummy_Regulatory_prob), 
        rep(1/11, 11),
        rep(1/208, 208),
        valueseq[i] * expert_hi
        )
      
      pred<- simbetas %*% xhyp_HI
      EVs_HI[i]<-quantile(pred,0.5)
      lower_HI[i]<-quantile(pred,0.025)
      upper_HI[i]<-quantile(pred,0.975)   } 

    out<-cbind(EVs_LO,lower_LO,upper_LO, EVs_HI, lower_HI, upper_HI)
    out

}

preds_A1<-predfunction_year(A1, A1_se) 
preds_A2<-predfunction_year(A2, A2_se) 
preds_A3<-predfunction_year(A3, A3_se) 
preds_A4<-predfunction_year(A4, A4_se) 
preds_A5<-predfunction_year(A5, A5_se) 



rm(A1)
rm(A2)
rm(A3)
rm(A4)
rm(A5)

rm(A1_se)
rm(A2_se)
rm(A3_se)
rm(A4_se)
rm(A5_se)






B1 <- lm(LOGtimeWD1 ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
B1_se<- cluster.vcov(B1, d2$DEPENDENCIA)

B2 <- lm(response_inexistencia ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
B2_se<- cluster.vcov(B2, d2$DEPENDENCIA)

B3 <- lm(R_trueinex ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
B3_se<- cluster.vcov(B3, d2$DEPENDENCIA)

B4 <- lm(R_lessthanhalf ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
B4_se<- cluster.vcov(B4, d2$DEPENDENCIA)

B5 <- lm(R_noncomply ~ 
  expert_combined*agency_expert_combined_mean_lag+
  logwordcount+readability +attach.inc+ medium_nonelec + inapproprequest + S4_num_info_queries_prob + S4_num_areas_prob + S5_distinct_reqs_related_prob + S7_dummy_Datum_prob + S7_dummy_Data_prob + S7_dummy_Database_prob + S7_dummy_Document_prob + S7_dummy_MultipleDocuments_prob + S8_dummy_Activities_prob + S8_dummy_Budget_prob + S8_dummy_Evaluation_prob + S8_dummy_ExternalContracts_prob + S8_dummy_InstStruc_prob + S8_dummy_Other_prob + S8_dummy_Regulatory_prob +
  as.factor(MONTH) + as.factor(DEPENDENCIA), data=d2)
B5_se<- cluster.vcov(B5, d2$DEPENDENCIA)




# confirm number of group fixed effects
length(names(coef(B1))[grep("DEPENDENCIA",names(coef(B1)))]) #[1] 205
length(names(coef(B2))[grep("DEPENDENCIA",names(coef(B2)))]) #205
length(names(coef(B4))[grep("DEPENDENCIA",names(coef(B4)))]) #205
length(names(coef(B5))[grep("DEPENDENCIA",names(coef(B5)))]) #205


predfunction_expertlag<-function(model, vcov){

    pe<-coef(model)
    vc<-vcov
    sims <- 1000
    simbetas <- mvrnorm(sims,pe,vc)

    valueseq<- seq(min(d2$agency_expert_combined_mean_lag, na.rm=T), max(d2$agency_expert_combined_mean_lag, na.rm=T), length.out=20)
    EVs_LO<-rep(NA,length(valueseq))
    lower_LO<-rep(NA,length(valueseq))
    upper_LO<-rep(NA,length(valueseq))
    EVs_HI<-rep(NA,length(valueseq))
    lower_HI<-rep(NA,length(valueseq))
    upper_HI<-rep(NA,length(valueseq))


    for (i in 1:length(valueseq)){
        xhyp_LO<-c( 1, 
        expert_combined=expert_lo,
        agency_expert_combined_mean_lag=valueseq[i],
        logwordcount = mean(d2$logwordcount),
        readability = mean(d2$readability),
        attach.inc = mean(d2$attach.inc),
        medium_nonelec = mean(d2$medium_nonelec),
        inapproprequest = mean(d2$inapproprequest),
        S4_num_info_queries_prob = mean(d2$S4_num_info_queries_prob),
        S4_num_areas_prob = mean(d2$S4_num_areas_prob),
        S5_distinct_reqs_related_prob = mean(d2$S5_distinct_reqs_related_prob),
        S7_dummy_Datum_prob = mean(d2$S7_dummy_Datum_prob),
        S7_dummy_Data_prob = mean(d2$S7_dummy_Data_prob),
        S7_dummy_Database_prob = mean(d2$S7_dummy_Database_prob),
        S7_dummy_Document_prob = mean(d2$S7_dummy_Document_prob),
        S7_dummy_MultipleDocuments_prob = mean(d2$S7_dummy_MultipleDocuments_prob),
        S8_dummy_Activities_prob = mean(d2$S8_dummy_Activities_prob),
        S8_dummy_Budget_prob = mean(d2$S8_dummy_Budget_prob),
        S8_dummy_Evaluation_prob = mean(d2$S8_dummy_Evaluation_prob),
        S8_dummy_ExternalContracts_prob = mean(d2$S8_dummy_ExternalContracts_prob),
        S8_dummy_InstStruc_prob = mean(d2$S8_dummy_InstStruc_prob),
        S8_dummy_Other_prob = mean(d2$S8_dummy_Other_prob),
        S8_dummy_Regulatory_prob = mean(d2$S8_dummy_Regulatory_prob),  
        rep(1/11, 11),
        rep(1/205, 205),
        valueseq[i] * expert_lo
        )
      
      pred<- simbetas %*% xhyp_LO
      EVs_LO[i]<-quantile(pred,0.5)
      lower_LO[i]<-quantile(pred,0.025)
      upper_LO[i]<-quantile(pred,0.975)   } 



    for (i in 1:length(valueseq)){
        xhyp_HI<-c( 1, 
        expert_combined=expert_hi,
        agency_expert_combined_mean_lag=valueseq[i],
        logwordcount = mean(d2$logwordcount),
        readability = mean(d2$readability),
        attach.inc = mean(d2$attach.inc),
        medium_nonelec = mean(d2$medium_nonelec),
        inapproprequest = mean(d2$inapproprequest),
        S4_num_info_queries_prob = mean(d2$S4_num_info_queries_prob),
        S4_num_areas_prob = mean(d2$S4_num_areas_prob),
        S5_distinct_reqs_related_prob = mean(d2$S5_distinct_reqs_related_prob),
        S7_dummy_Datum_prob = mean(d2$S7_dummy_Datum_prob),
        S7_dummy_Data_prob = mean(d2$S7_dummy_Data_prob),
        S7_dummy_Database_prob = mean(d2$S7_dummy_Database_prob),
        S7_dummy_Document_prob = mean(d2$S7_dummy_Document_prob),
        S7_dummy_MultipleDocuments_prob = mean(d2$S7_dummy_MultipleDocuments_prob),
        S8_dummy_Activities_prob = mean(d2$S8_dummy_Activities_prob),
        S8_dummy_Budget_prob = mean(d2$S8_dummy_Budget_prob),
        S8_dummy_Evaluation_prob = mean(d2$S8_dummy_Evaluation_prob),
        S8_dummy_ExternalContracts_prob = mean(d2$S8_dummy_ExternalContracts_prob),
        S8_dummy_InstStruc_prob = mean(d2$S8_dummy_InstStruc_prob),
        S8_dummy_Other_prob = mean(d2$S8_dummy_Other_prob),
        S8_dummy_Regulatory_prob = mean(d2$S8_dummy_Regulatory_prob), 
        rep(1/11, 11),
        rep(1/205, 205),
        valueseq[i] * expert_hi
        )
      
      pred<- simbetas %*% xhyp_HI
      EVs_HI[i]<-quantile(pred,0.5)
      lower_HI[i]<-quantile(pred,0.025)
      upper_HI[i]<-quantile(pred,0.975)   } 


    out<-cbind(EVs_LO,lower_LO,upper_LO, EVs_HI, lower_HI, upper_HI)
    out

}

preds_B1<-predfunction_expertlag(B1, B1_se) 
preds_B2<-predfunction_expertlag(B2, B2_se) 
preds_B3<-predfunction_expertlag(B3, B3_se) 
preds_B4<-predfunction_expertlag(B4, B4_se) 
preds_B5<-predfunction_expertlag(B5, B5_se) 


rm(B1)
rm(B2)
rm(B3)
rm(B4)
rm(B5)

rm(B1_se)
rm(B2_se)
rm(B3_se)
rm(B4_se)
rm(B5_se)





#Figure 2: Visualizations of key results from Models 1 and 3 in Table 2 and Models 1 and 3 in Table 3.

valueseq<- seq(min(d2$agency_expert_combined_mean_lag, na.rm=T), max(d2$agency_expert_combined_mean_lag, na.rm=T), length.out=20)

pdf(file="Figure2.pdf", height= 9, width=14 )

par(mfrow=c(2,2))

plot(2003:2019, exp(preds_A1[,1]), type="l", lwd=3, col="gray70", ylim=range(exp(preds_A1)), xaxt="n",
  xlab="Year", ylab="Predicted Working Days to Response", main="Differences by Request Expertise in Working Days to Response\nVarying Over Time")
lines(2003:2019, exp(preds_A1[,2]), lwd=1, lty=1, col="gray70")
lines(2003:2019, exp(preds_A1[,3]), lwd=1, lty=1, col="gray70")
lines(2003:2019, exp(preds_A1[,4]), lwd=3, col="black")
lines(2003:2019, exp(preds_A1[,5]), lwd=1, lty=1, col="black")
lines(2003:2019, exp(preds_A1[,6]), lwd=1, lty=1, col="black")
axis(side=1, at=2003:2019)
legend(x=2010, y=15, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))



plot(2003:2019, preds_A3[,1], type="l", lwd=3, col="gray70", ylim=range(preds_A3), xaxt="n",
  xlab="Year", ylab="Predicted Probability of True 'Inexistencia'", main="Differences by Request Expertise in True 'Inexistencia' Responses\nVarying Over Time")
lines(2003:2019, preds_A3[,2], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A3[,3], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A3[,4], lwd=3, col="black")
lines(2003:2019, preds_A3[,5], lwd=1, lty=1, col="black")
lines(2003:2019, preds_A3[,6], lwd=1, lty=1, col="black")
axis(side=1, at=2003:2019)
legend(x=2010, y=0.1, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(valueseq, exp(preds_B1[,1]), type="l", lwd=3, col="gray70", ylim=range(exp(preds_B1)), 
  xlab="Agency's Previous Year Average Expert Request Score", ylab="Predicted Working Days to Response", main="Differences by Request Expertise in Working Days to Response\nVarying Over Agency's Previous Year Average Expert Request Score")
lines(valueseq, exp(preds_B1[,2]), lwd=1, lty=1, col="gray70")
lines(valueseq, exp(preds_B1[,3]), lwd=1, lty=1, col="gray70")
lines(valueseq, exp(preds_B1[,4]), lwd=3, col="black")
lines(valueseq, exp(preds_B1[,5]), lwd=1, lty=1, col="black")
lines(valueseq, exp(preds_B1[,6]), lwd=1, lty=1, col="black")
legend(x=0.1, y=18, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(valueseq, preds_B3[,1], type="l", lwd=3, col="gray70", ylim=range(preds_B3), 
  xlab="Agency's Previous Year Average Expert Request Score", ylab="Predicted Probability of True 'Inexistencia'", main="Differences by Request Expertise in True 'Inexistencia' Responses\nVarying Over Agency's Previous Year Average Expert Request Score")
lines(valueseq, preds_B3[,2], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B3[,3], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B3[,4], lwd=3, col="black")
lines(valueseq, preds_B3[,5], lwd=1, lty=1, col="black")
lines(valueseq, preds_B3[,6], lwd=1, lty=1, col="black")
legend(x=0.1, y=0.35, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))

dev.off()





#Figure G.4: Visualizations of key results from Models 2, 4, and 5 in Table 2 and Models 2, 4, and 5 in Table 3. 

pdf(file="FigureG4.pdf", height= 9, width=16 )

par(mfrow=c(2,3))

plot(2003:2019, preds_A2[,1], type="l", lwd=3, col="gray70", ylim=range(preds_A2), xaxt="n",
  xlab="Year", ylab="Predicted Probability of Official 'Inexistencia'", main="Differences by Request Expertise in Official 'Inexistencia' Responses\nVarying Over Time")
lines(2003:2019, preds_A2[,2], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A2[,3], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A2[,4], lwd=3, col="black")
lines(2003:2019, preds_A2[,5], lwd=1, lty=1, col="black")
lines(2003:2019, preds_A2[,6], lwd=1, lty=1, col="black")
axis(side=1, at=2003:2019)
legend(x=2010, y=0.1, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))



plot(2003:2019, preds_A4[,1], type="l", lwd=3, col="gray70", ylim=range(preds_A4), xaxt="n",
  xlab="Year", ylab="Predicted Probability of Uninformative Response", main="Differences by Request Expertise in Uninformative Responses\nVarying Over Time")
lines(2003:2019, preds_A4[,2], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A4[,3], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A4[,4], lwd=3, col="black")
lines(2003:2019, preds_A4[,5], lwd=1, lty=1, col="black")
lines(2003:2019, preds_A4[,6], lwd=1, lty=1, col="black")
axis(side=1, at=2003:2019)
legend(x=2010, y=0.63, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(2003:2019, preds_A5[,1], type="l", lwd=3, col="gray70", ylim=range(preds_A5), xaxt="n",
  xlab="Year", ylab="Predicted Probability of Noncompliant Response", main="Differences by Request Expertise in Noncompliant Responses\nVarying Over Time")
lines(2003:2019, preds_A5[,2], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A5[,3], lwd=1, lty=1, col="gray70")
lines(2003:2019, preds_A5[,4], lwd=3, col="black")
lines(2003:2019, preds_A5[,5], lwd=1, lty=1, col="black")
lines(2003:2019, preds_A5[,6], lwd=1, lty=1, col="black")
axis(side=1, at=2003:2019)
legend(x=2010, y=0.1, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(valueseq, preds_B2[,1], type="l", lwd=3, col="gray70", ylim=range(preds_B2), 
  xlab="Agency's Previous Year Average Expert Request Score", ylab="Predicted Probability of Official 'Inexistencia'", main="Differences by Request Expertise in Official 'Inexistencia' Responses\nVarying Over Agency's Previous Year Average Expert Request Score")
lines(valueseq, preds_B2[,2], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B2[,3], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B2[,4], lwd=3, col="black")
lines(valueseq, preds_B2[,5], lwd=1, lty=1, col="black")
lines(valueseq, preds_B2[,6], lwd=1, lty=1, col="black")
legend(x=0.1, y=-0.05, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(valueseq, preds_B4[,1], type="l", lwd=3, col="gray70", ylim=range(preds_B4), 
  xlab="Agency's Previous Year Average Expert Request Score", ylab="Predicted Probability of Uninformative Response", main="Differences by Request Expertise in Uninformative Responses\nVarying Over Agency's Previous Year Average Expert Request Score")
lines(valueseq, preds_B4[,2], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B4[,3], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B4[,4], lwd=3, col="black")
lines(valueseq, preds_B4[,5], lwd=1, lty=1, col="black")
lines(valueseq, preds_B4[,6], lwd=1, lty=1, col="black")
legend(x=0.1, y=0.3, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))


plot(valueseq, preds_B5[,1], type="l", lwd=3, col="gray70", ylim=range(preds_B5), 
  xlab="Agency's Previous Year Average Expert Request Score", ylab="Predicted Probability of Noncompliant Response", main="Differences by Request Expertise in Noncompliant Responses\nVarying Over Agency's Previous Year Average Expert Request Score")
lines(valueseq, preds_B5[,2], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B5[,3], lwd=1, lty=1, col="gray70")
lines(valueseq, preds_B5[,4], lwd=3, col="black")
lines(valueseq, preds_B5[,5], lwd=1, lty=1, col="black")
lines(valueseq, preds_B5[,6], lwd=1, lty=1, col="black")
legend(x=0.1, y=0.3, lty=c(1,1), lwd=c(2,2), col=c("black","gray70"), bty="n", cex=0.9,
  legend=c("Expert Request Score = 1 SD above Mean","Expert Request Score = 1 SD below Mean"))

dev.off()







### Nonlinear Interaction Test Figures

#Figures for Figure H.1 and Figure H.2

out1 <- interflex(Y = "LOGtimeWD1", D = "expert_combined", X = "year_ind", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out1, Ylabel="Log Days-to-Response", Dlabel="Expert Request Score", Xlabel="Years", file="interflex_plot_1a.pdf")
print(out1$tests$p.wald)
rm(out1)

out2 <- interflex(Y = "response_inexistencia", D = "expert_combined", X = "year_ind", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out2, Ylabel="Official Inexistencia", Dlabel="Expert Request Score", Xlabel="Years", file="interflex_plot_1b.pdf")
print(out2$tests$p.wald)
rm(out2)

out3 <- interflex(Y = "R_trueinex", D = "expert_combined", X = "year_ind", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out3, Ylabel="True Inexistencia", Dlabel="Expert Request Score", Xlabel="Years", file="interflex_plot_1c.pdf")
print(out3$tests$p.wald)
rm(out3)

out4 <- interflex(Y = "R_lessthanhalf", D = "expert_combined", X = "year_ind", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out4, Ylabel="Uninformative", Dlabel="Expert Request Score", Xlabel="Years", file="interflex_plot_1d.pdf")
print(out4$tests$p.wald)
rm(out4)

out5 <- interflex(Y = "R_noncomply", D = "expert_combined", X = "year_ind", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out5, Ylabel="Noncompliant", Dlabel="Expert Request Score", Xlabel="Years", file="interflex_plot_1e.pdf")
print(out5$tests$p.wald)
rm(out5)



out1a <- interflex(Y = "LOGtimeWD1", D = "expert_combined", X = "agency_expert_combined_mean_lag", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out1a, Ylabel="Log Days-to-Response", Dlabel="Expert Request Score", Xlabel="Agency Lagged Avg. Expertise", file="interflex_plot_2a.pdf")
print(out1a$tests$p.wald)
rm(out1a)

out2a <- interflex(Y = "response_inexistencia", D = "expert_combined", X = "agency_expert_combined_mean_lag", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out2a, Ylabel="Official Inexistencia", Dlabel="Expert Request Score", Xlabel="Agency Lagged Avg. Expertise", file="interflex_plot_2b.pdf")
print(out2a$tests$p.wald)
rm(out2a)

out3a <- interflex(Y = "R_trueinex", D = "expert_combined", X = "agency_expert_combined_mean_lag", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out3a, Ylabel="True Inexistencia", Dlabel="Expert Request Score", Xlabel="Agency Lagged Avg. Expertise", file="interflex_plot_2c.pdf")
print(out3a$tests$p.wald)
rm(out3a)

out4a <- interflex(Y = "R_lessthanhalf", D = "expert_combined", X = "agency_expert_combined_mean_lag", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out4a, Ylabel="Uninformative", Dlabel="Expert Request Score", Xlabel="Agency Lagged Avg. Expertise", file="interflex_plot_2d.pdf")
print(out4a$tests$p.wald)
rm(out4a)

out5a <- interflex(Y = "R_noncomply", D = "expert_combined", X = "agency_expert_combined_mean_lag", 
  Z = c("logwordcount","readability","attach.inc","medium_nonelec","inapproprequest","S4_num_info_queries_prob","S4_num_areas_prob","S5_distinct_reqs_related_prob","S7_dummy_Datum_prob","S7_dummy_Data_prob","S7_dummy_Database_prob","S7_dummy_Document_prob","S7_dummy_MultipleDocuments_prob","S8_dummy_Activities_prob","S8_dummy_Budget_prob","S8_dummy_Evaluation_prob","S8_dummy_ExternalContracts_prob","S8_dummy_InstStruc_prob","S8_dummy_Other_prob","S8_dummy_Regulatory_prob"), 
  FE = c("MONTH","DEPENDENCIA"), cl = "DEPENDENCIA", 
  data = d2, na.rm=T,
  estimator = "binning", vcov.type = "cluster", main = "Marginal Effects", ylim = c(-15, 15))

plot(out5a, Ylabel="Noncompliant", Dlabel="Expert Request Score", Xlabel="Agency Lagged Avg. Expertise", file="interflex_plot_2e.pdf")
print(out5a$tests$p.wald)
rm(out5a)








### Analyses of Demand-Side Feedback Effects

d5<-d1[d1$RESPUESTA!="Sin respuesta" & !is.na(d1$LOGtimeWD1) & !is.na(d1$CODE_COMBINED) & d1$YEAR<2017,]

# create a muni x year dataframe
um<-unique(d5$CODE_COMBINED)
muni_df<- as.data.frame(cbind(
  rep(um, each=length(2003:2016)),
  rep(2003:2016, length(um))
  ))
names(muni_df)<-c("CODE_COMBINED","YEAR")
muni_df$YEAR<-as.numeric(muni_df$YEAR)

#avg for each response variable
muni_df$requestvolume<-NA
muni_df$LOGtimeWD1_MUNImean<-NA
muni_df$response_inexistencia_MUNImean<-NA
muni_df$R_trueinex_MUNImean<-NA
muni_df$R_lessthanhalf_MUNImean<-NA
muni_df$R_noncomply_MUNImean<-NA

timestamp()
for(i in 1:nrow(muni_df)){
  tempdata<-d5[d5$CODE_COMBINED==muni_df$CODE_COMBINED[i] & d5$YEAR==muni_df$YEAR[i],]
  muni_df$requestvolume[i]<-nrow(tempdata)
  muni_df$LOGtimeWD1_MUNImean[i]<-mean(tempdata$LOGtimeWD1, na.rm=T)
  muni_df$response_inexistencia_MUNImean[i]<-mean(tempdata$response_inexistencia, na.rm=T)
  muni_df$R_trueinex_MUNImean[i]<-mean(tempdata$R_trueinex, na.rm=T)
  muni_df$R_lessthanhalf_MUNImean[i]<-mean(tempdata$R_lessthanhalf, na.rm=T)
  muni_df$R_noncomply_MUNImean[i]<-mean(tempdata$R_noncomply, na.rm=T)

  if(i %in% seq(100,30000,1000)) print(i)
}
timestamp()


#lagging variables
muni_df$requestvolume_lag<-lagger(muni_df$requestvolume, muni_df$CODE_COMBINED, muni_df$YEAR, 1)
muni_df$LOGtimeWD1_MUNImean_lag<-lagger(muni_df$LOGtimeWD1_MUNImean, muni_df$CODE_COMBINED, muni_df$YEAR, 1)
muni_df$response_inexistencia_MUNImean_lag<-lagger(muni_df$response_inexistencia_MUNImean, muni_df$CODE_COMBINED, muni_df$YEAR, 1)
muni_df$R_trueinex_MUNImean_lag<-lagger(muni_df$R_trueinex_MUNImean, muni_df$CODE_COMBINED, muni_df$YEAR, 1)
muni_df$R_lessthanhalf_MUNImean_lag<-lagger(muni_df$R_lessthanhalf_MUNImean, muni_df$CODE_COMBINED, muni_df$YEAR, 1)
muni_df$R_noncomply_MUNImean_lag<-lagger(muni_df$R_noncomply_MUNImean, muni_df$CODE_COMBINED, muni_df$YEAR, 1)


# merging these back into main data
# using a temporary smaller dataframe without extraneous variables so that this will run faster, then merging back in
d1_temp<-d1[d1$YEAR %in% 2004:2016, c("FOLIO","YEAR","CODE_COMBINED")]

d1_temp$requestvolume_lag<-NA
d1_temp$LOGtimeWD1_MUNImean_lag<-NA
d1_temp$response_inexistencia_MUNImean_lag<-NA
d1_temp$R_trueinex_MUNImean_lag<-NA
d1_temp$R_lessthanhalf_MUNImean_lag<-NA
d1_temp$R_noncomply_MUNImean_lag<-NA

# skip rows of the municipality-year dataframe where request volume is zero, to speed merge
muni_df_temp<-muni_df[muni_df$requestvolume_lag>0 & !is.na(muni_df$requestvolume_lag),]
dim(muni_df_temp)

# this can take up to an hour
timestamp()
for(i in 1:nrow(muni_df_temp)){
  yeartemp<-muni_df_temp$YEAR[i]
  munitemp<-muni_df_temp$CODE_COMBINED[i]

  d1_temp$requestvolume_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$requestvolume_lag[i]
  d1_temp$LOGtimeWD1_MUNImean_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$LOGtimeWD1_MUNImean_lag[i]
  d1_temp$response_inexistencia_MUNImean_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$response_inexistencia_MUNImean_lag[i]
  d1_temp$R_trueinex_MUNImean_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$R_trueinex_MUNImean_lag[i]
  d1_temp$R_lessthanhalf_MUNImean_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$R_lessthanhalf_MUNImean_lag[i]
  d1_temp$R_noncomply_MUNImean_lag[d1_temp$YEAR==yeartemp & d1_temp$CODE_COMBINED==munitemp]<- muni_df_temp$R_noncomply_MUNImean_lag[i]

  if(i %in% seq(100,13000,100)) print(i)

}
timestamp()


# merging back together

d1_temp$CODE_COMBINED<-NULL
d1_temp$YEAR<-NULL

d1_new<-merge(d1, d1_temp, by="FOLIO")
d2_new<-d1_new[d1_new$YEAR<2017, ] 



#Table F.3: Linear models of request expertise, for requests 2003-2016. 

m1<-felm(expert_combined ~ LOGtimeWD1_MUNImean_lag  | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)
m2<-felm(expert_combined ~ LOGtimeWD1_MUNImean_lag  | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)

m3<-felm(expert_combined ~ response_inexistencia_MUNImean_lag  | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)
m4<-felm(expert_combined ~ response_inexistencia_MUNImean_lag  | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)

m5<-felm(expert_combined ~ R_trueinex_MUNImean_lag  | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)
m6<-felm(expert_combined ~ R_trueinex_MUNImean_lag  | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)

m7<-felm(expert_combined ~ R_lessthanhalf_MUNImean_lag  | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)
m8<-felm(expert_combined ~ R_lessthanhalf_MUNImean_lag  | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)

m9<-felm(expert_combined ~ R_noncomply_MUNImean_lag  | 0 | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)
m10<-felm(expert_combined ~ R_noncomply_MUNImean_lag  | CODE_COMBINED | 0 | CODE_COMBINED, exactDOF=TRUE,psdef=FALSE, data=d2_new)

texreg(list(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10), digits=4,  stars = c(0.01, 0.05, 0.1), omit.coef=c("as.factor"))






















